clc;clear;clearvars;
% 随机生成5个数据
num_initial = 5;
num_vari = 5;
% 搜索区间
upper_bound = 32;
lower_bound = -32;
iter = 3000;
% 随机生成5个数据,并获得其评估值
sample_x = lhsdesign(num_initial, num_vari).*(upper_bound - lower_bound) + lower_bound.*ones(num_initial, num_vari);
sample_y = F2(sample_x);
Fmin = zeros(iter, 1);
aver_Fmin = zeros(iter, 1);
num1 = 0;

for n = 1 : 100
    k = 1;
    % 初始化一些参数
    pbestx = sample_x;
    pbesty = sample_y;
    [fmin, gbest] = min(pbesty);
    fprintf("iter 0 fmin: %.4f\n", fmin);
    for i = 1 : iter
        % pso更新下一步的位置,这里可以设置一下超过搜索范围的就设置为边界
        r = rand;
        delta = zeros(num_initial, num_vari);
        for i1 = 1 : num_initial
            for i2 = 1 : num_vari
                k1 = randi(num_initial, 1, 1);
                k2 = randi(num_initial, 1, 1);
                delta(i1,i2) = 1 * rand * abs(pbestx(k1,i2) - pbestx(k2,i2)) * (exp(pbesty(i1,1) - pbesty(gbest,1)));
            end
        end
        if r > 0.5
            x = pbestx;
        else
            x = normrnd((pbestx + pbestx(gbest, :)) ./ 2,abs(pbestx - pbestx(gbest, :)) + delta);
        end

        % The adaptive mutation
        pc = 0.2 / (1 + 1.2^(20 - num1));
        for i4 = 1 : num_initial
            r3 = rand;
            if r3 < pc
                u_n = ceil(pc * num_vari);
                for i5 = 1 : u_n
                    l = randi(num_vari, 1, 1);
                    x(i4, l) = rand * (upper_bound - lower_bound) + lower_bound;
                end
            end
        end

        x(x > upper_bound) = upper_bound;
        x(x < lower_bound) = lower_bound;
        y = F2(x);
        % 更新每个单独个体最佳位置
        pbestx(y < pbesty, :) = x(y < pbesty, :);
        pbesty = F2(pbestx);
        % 更新所有个体最佳位置
        [fmin1, gbest] = min(pbesty);
        if fmin1 < fmin
            num1 = 0;
        else
            num1 = num1 + 1;
        end
        fmin = fmin1;
        if mod(i,100) == 0
            fprintf("iter %d fmin: %.4f\n", i, fmin);
        end
        Fmin(k, 1) = fmin;
        k = k +1;
    end
    aver_Fmin = aver_Fmin + Fmin;
end
aver_Fmin = aver_Fmin ./ 100;
disp(pbestx(gbest, :));
plot(aver_Fmin);